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Cosmological perturbations in the DGP braneworld: numeric solution 
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We solve for the behaviour of cosmological perturbations in the Dvali-Gabadadze-Porrati (DGP) 
braneworld model using a new numerical method. Unlike some other approaches in the literature, 
our method uses no approximations other than linear theory and is valid on large scales. We exam- 
ine the behaviour of late-universe density perturbations for both the self-accelerating and normal 
branches of DGP cosmology. Our numerical results can form the basis of a detailed comparison 
between the DGP model and cosmological observations. 



I. INTRODUCTION 

In the braneworld paradigm our universe is a lower- 
dimensional object embedded in a higher-dimensional 
bulk spacetime. Standard Model fields are assumed to 
be confined to the brane while gravity is allowed to prop- 
agate in the bulk. The Dvali-Gabadadze-Porrati (DGP) 
[J Q model postulates that we live in a 4-dimensional 
hypersurface in a 5-dimensional Minkowski bulk. Gen- 
eral Relativity (GR) is recovered at small scales (smaller 
than the crossover scale Tc) due to the inclusion of an 
induced gravity term in the action. 

It was quickly realized that this model has two dis- 
tinct classes of cosmological solutions Q. One of them 
exhibits accelerated expansion at late times without the 
need to include any exotic cosmological fluids, such as 
dark energy, or any brane tension that acts as an effective 
4-dimensional cosmological constant. Hence, this branch 
of solutions is called "self-accelerating" . Several attempts 
to confront the self-accelerating universe with observa- 
tions have been made [1, H, @ (also see 0, review] and 
references therein). To explain the observed acceleration 
we require Tc ~ , where Hq is the current value of 
the Hubble parameter. It is expected that structure for- 
mation will help to distinguish the self-accelerating DGP 
universe from dark energy models based on 4-dimensional 
GR. This is because the growth of cosmological pertur- 
bations is very sensitive to the existence of an extra di- 
mension. A full 5-dimensional treatment is required to 
model these perturbations, which is why obtaining ob- 
servational predictions for the behaviour of fluctuations 
in the DGP model is technically challenging. 

Several authors have considered the problem of the dy- 
namics of perturbations in the DGP model, but they have 
all relied on some sort of approximation or simplifying 
ansatz. Two examples of this are the quasi-static (QS) 
approximation developed by Koyama and Maartens [3| 
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and the dynamical scaling (DS) ansatz proposed by Saw- 
icki et al. Q. The former approximation scheme solves 
the perturbative equations of motion by focussing on the 
extreme subhorizon regime. In contrast the DS method, 
which assumes that perturbations evolve as power laws 
of the scale factor with time-varying power law indices, 
is supposed to be valid on all scales. It has been shown 
that the DS solution approaches the QS solution for sub- 
horizon perturbations. 

In this paper we present a complete numerical analy- 
sis of the evolution of scalar perturbations in the DGP 
model. Mathematically, the problem involves the solu- 
tion of a partial differential equation in the bulk cou- 
pled to an ordinary differential equation on the brane. 
A numerical method for dealing with such systems has 
previously been developed for cosmological perturbations 
in the Randall-Sundrum (RS) model [iOilllI- However, 
the DGP problem is more complicated than the RS case 
due to a non-local boundary condition on the bulk field. 
Hence, the algorithm used in this paper represents a sig- 
nificant generalization of the one used in Refs. [Tol. IT]|. 

Unfortunately, some theoretical issues cast doubt 
on the validity of the self-accelerating DGP solutions. 
Specifically, the existence of a perturbative ghost perhaps 
suggest s that this solution cannot describe our Universe 
Il3| , though there has been a debate on the physical 
implication of the ghost mode (see (T^ . review] and ref- 
erences therein). However, as alluded to above there is 
another "normal" branch of solutions in the DGP model. 
We cannot explain the late time accelerated expansion of 
the Universe using the normal branch without including 
an effective cosmological constant induced by the brane 
tension a. However, by allowing a non-zero cr, normal 
branch solutions can mimic dark energy models with the 
equation of state w smaller than —1 [l^, [l^. Cosmo- 
logical constraints on the background dynamics of these 
models have been studied [13, [Sj ■ Unlike 4-dimensional 
models that realize w < —1 by the introduction of a 
phantom field, the normal branch of DGP cosmology is 
ghost-free [l^. This unique feature is what motivates 
us to numerically study the perturbations of the normal 
branch of DGP cosmology in the penultimate section of 
this paper. During the preparation of this manuscript. 
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we learned that the behaviour of perturbations in the 
normal branch has been independently obtained by Song 
pol l using the DS method. 

The structure of the paper is as follows: The back- 
ground cosmology of the DGP model is discussed in SjlTl 
In we express the equations of motion for scalar and 
tensor perturbations of the DGP model in the dimension- 
less canonical form introduced in §111 Al The numerical 
method used to solve these canonical equations of motion 
is developed in ^IVl Our algorithm is tested in ^Vl where 
we numerically recover analytic results for the behaviour 
of tensor perturbations in matter-free DGP models. In 
i jVIl and Will we solve the scalar perturbations problem 
in the self-accelerating and normal branches, respectively. 
Finally, Willi is reserved for our conclusions. 



II. BACKGROUND SOLUTION 

A. Field equations and junction conditions 

We consider a 5-dimensional manifold A4 with metric 
gab, and covered by coordinates {X^j^^g. The manifold 
has a 4-dimensional brane boundary dMh with intrinsic 
metric jap, and covered by coordinates {x'^}a=o- The 
action of the model is 

a). (1) 




Here, a is the brane tension and Cm is the matter La- 
grangian. We impose the Z2 symmetry that the bulk is 
mirror symmetric about the brane. ^ The field equations 
satisfied in A4 are simply 



(2) 



We write the brane normal pointing into M. as n°. We 
find that the brane's extrinsic curvature, 



(3) 



must satisfy 



Kaf3 - Kja0 - rcG^^l = -^Kl{Taf3 - Crffa/s), (4) 



where we have defined the cross-over distance Tc by 



Tr. = 



(5) 



We assume that the stress-energy tensor of the brane 
matter. 



T, 



Q/3 



^7 



Q/3 



{p+p)UaUp+pgaP, (6) 



is of the perfect fluid form. Note that as in GR, the stress 
energy tensor is conserved, V^Tap = 0. Using the Gauss- 
Codazzi equations, it is possible to re-write the junction 
conditions ^ as the effective Einstein equations 



^(4) 



(7) 



where 



n, 



1 rri rri Q 

" 4 ^ fia-^ V 



—TT 

12^ 



T„ 



^9af}- Ki'^G^^}, (8) 



and £^1, is the trace-free projection of the 5-dimensional 
Weyl tensor. 

An important tool that is often used to analyze 
braneworld models is the "Gaussian-normal" coordi- 
nates. These are constructed by looking at the spatial 
geodesies that extend perpendicularly from dMh into A^. 
Gaussian-normal coordinates are then given by 
where y is the afHne parameter along these geodesies and 
are 4-dimensional coordinates on the family of hyper- 
surfaces tangent to the brane. Without loss of generality, 
we can set the brane to be at ?/ = 0. Most importantly, 
the derivative of any bulk scalar quantity with respect to 
y on the brane corresponds to the normal derivative: 



dy 



(9) 



Below, (dy • • • )b will always be understood to be the nor- 
mal derivative of some quantity evaluated at the brane. 



B. Bulk geometry and brane trajectory 

One solution of the above field equations makes use of 
the following 5-dimensional flat metric with Rabcd = 0: 



ds^ = gabdX'^dX'' = -r^ du dv + dy? 



(10) 



Here, u and v are dimensionless null coordinates. The 
brane is defined parametrically by (u, v) — {u\j(t), Wb(t)), 
where 



v^(t) = a{t), Ub{t) = \ [ 

ri Jo a.[x) 



(11) 



Here, a{t) is the scale factor of the brane universe, nor- 
malized to unity today, and t is the proper time along 
the brane. The latter implies that the following induced 
line element on the brane is of the FRW form: 



^ Technically, M refers to one half of the total bulk spacetime. 



dsl = -fafidx^'dx'^ = -dt^ + a^{t) dx^ (12) 
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We have selected the u coordinate such that 7Xb(0) = 
and the brane moves in the direction of increasing u and 
V. The 5-velocity of the brane and brane normal 
are given by 



dt = u-'da 



dx" d 



dy = n'^da = e 
These satisfy 



dt dx" 
1 



-rrdu 



ddv, 



ria 



u°-Ua = -1, u°-na 0, n'^ria = 1. 



(13a) 
(13b) 

(14) 



The e = ±1 parameter in the definition of the brane 
normal (|13b[) reflects the fact that when we impose the 
Z2 symmetry across the brane, we have two choices for 
the half of the bulk manifold M we discard. 

The junction conditions (H]) can be used to determine 
the brane dynamics, which are governed by 



H 



a 
a 



1 

2r7 



e+\ll + -Klrl{p + a) 



w)pH, 



dH rcKl(l + w)pH 



(15a) 

(15b) 
(15c) 



Note that when e = +1, we have that Hvc ~ 1 when the 
density of brane matter is small |p + (t| ^ n^^r^^. This 
implies a late-time accelerating universe, which is why 
the e = +1 case is called the self- accelerating branch and 
the e = — 1 case is called the normal branch. 



III. MASTER EQUATIONS GOVERNING 
PERTURBATIONS 



A. Dimensionless coordinates and canonical wave 
equations 

In this paper, we will consider the perturbations of the 
DGP model governed by a field defined on the bulk 
spacetime M coupled to a dynamical field A residing on 
the brane dMh- It is useful to decompose these fields 
into Fourier modes as follows: 



V'(m,u,x) = 
A(i,x) = 



d^k 
(27r)3/2 

d^k 
(27r)3/2 



Ak(i)e*'^-''. 



(16a) 
(16b) 



As usual for linear theory, the individual k modes are de- 
coupled from one another. In what follows, we will omit 
the k subscript from and Ak. That is, ip and A re- 
fer to the Fourier amplitudes of modes with wavevector 
k. We will also assume that a normalization has been 



selected such that the Fourier amplitudes are dimension- 
less. 

For a given mode, we define the "*" epoch as the mo- 
ment when a given mode crosses the Hubble horizon: 



k — H^ci^^ 



kk. 



(17) 



Then, we can define a set of normalized variables (deco- 
rated with hats): 



H = Hr,, 
2 2 

d — a/a-j,, 
Explicitly, wc have 



"-2 



k = 


krc/a^,, 


P = 


2 2 
K^r^p, 


i = 




y = 


y/rc, 


u — 




V = 


v/a*. 



1+ 3(/5 + <^) 



(18) 

(19a) 
(19b) 



It is also useful to define the 2-dimensional Cartesian 
coordinates 



n V + u 

X = T = — - — , 



(20) 



In terms of these coordinates, the dimensionless tangen- 
tial dt = U^Oa and normal dy = N^Sa derivatives to 
the brane are 



dt 



4 = 1 



Hd 

Hd^ 



Hd) 
J_\ d_ 
Hd) dr 



Hd) dz 



-T-+\Hd 



1 



d 



Hd) dz^ 



where 



1 ^ UaU"^, 1 = NaN"^, = UaN^, 



(21a) 
(21b) 

(22) 



and capital roman indices A, B — 0, 1 are raised and 
lowered with the flat 2-metric t^ab — diag(— 1, 1). 

The brane trajectory in the (r, z) plane is given by the 
solution of 



Hd. 



dd 
di 

dTk I (fj. ^ 1 
— - = — Ha + — — 
dt 2\ Hd 

dz^ 1 ( ^ 1 

— — — — Ha — 

dt 2\ Hd 



subject to the initial conditions 

a(0) = di, Tb(0) = \di, Zb(0) 



(23a) 
(23b) 
(23c) 

(24) 



If e = -|-1, the normal points in the direction of increas- 
ing z and the "bulk" corresponds to the the portion of 
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the (r, z) plane to the right of the brane; if e = — 1 the 
opposite is true. 

We win find below that the equations governing the 
bulk ifj = iIj{t, z) and brane field A — A(f) are of the 
form: 

0^{dl- dl + (25a) 
(9yV)b = AiA + AaS + Agi^b + Ui^'^, + Ag^, (25b) 
S' = AeA + A7S + Asi/^b + AgV-b + AioV'b, (25c) 
A' = S, (25d) 

where we have introduced the auxiliary field S, which 
corresponds to the time derivative of A. In (j25p . a prime 
' denotes the derivative with respect to the dimensionless 
proper time dj dt and i/'b — "iphit) — 4^(jh{t), Zh(t)) is the 
value of the bulk field on the brane. Also, V — V{t, z) 
is the bulk potential and the coefficients A„ = A„(a) are 
functions of the brane scale factor. All quantities ap- 
pearing in are dimensionless. We refer to (pS)) as 
the canonical form of the perturbative equations of mo- 
tion. 



B. Density perturbations in the late universe 

In this subsection we describe the formulae governing 
scalar perturbations in the late-time matter dominated 
universe. We take the matter content of the brane to be 
a dust fluid w = (i.e., cold dark matter): 



p (X a 



(26) 



In the 5-dimensional longitudinal gauge, scalar-type 
perturbations^ of the bulk geometry pUj) can be written 
as 

ds^ = — {du dv + Fuudu^ + 2Fuvdu dv + F^ydv^) 
+ rdfmdu + fvidv)dx' + v^{l + 2n)6,jdx'dx^ . (27) 

It can be shown that the dynamics of all of the per- 
turbative quantities in this expression can be derived 
from a single scalar bulk degree of freedom [2l|. After 
Fourier decomposition, we find that the mode amplitude 
— n{u,v) of this master field obeys 



= 



dudv 2v du Av"^ 



(28) 



We parameterize the fluctuations of the brane geometry 
by the two metric potentials $ and ^: 

dsl = -(1 + 21')rfi2 + a^{i + 2^)S,jdx'dx'. (29) 



^ That is, perturbations that can be derived from scalar potentials 
defined on the 3-dimension spatial sections. 



Finally, we write the perturbation of the brane fluid stress 
energy tensor as 

ST°o^-Sp, 6T°i^adiSq, ST'j^SpS'j. (30) 

One can construct the comoving density perturbation 
from these quantities as follows: 



pA ~ Sp — 3Ha 5q. 



(31) 



Both A and Q are gauge invariant quantities under a co- 
ordinate transformation on the brane and in the bulk. 
Note that A is the density contrast of the cold dark mat- 
ter only, not the the total density contrast. 

As shown in [23 | and demonstrated in Appendix [Aj 
one can use the effective Einstein equations on the brane 
and the detailed relationship between Q and the metric 
perturbations to flnd the following boundary condition 

for n-. 



3(e73fc2 -t- j^H'^a 



'"''»b + '^^-:y^^ A; (32) 



4iJa2 " 2fc2 

the following equation of motion for A: 

A + 2HA- ^^2^72 A = -'-^n^ 
I 4a 

and the following expressions for $ and ^: 

* = +^A ^il, - ^(^^+3g^f)7i^^^ 



(33) 



2fc2 4arc 



4a 



4rca3 



rib, 



(34) 



where fib = ^h{t) = ^{ub{t),Vh(t)). In these expres- 
sions, the dimensionless 7-factors are: 



71 



72 



73 



74 



2eHr^ 
2eHrc - 1 ' 

2erc{H - + 2eH^r^) 
H{2eHrc - 1)^ ' 

4erc(2ercff - 3iJ + 6ei/Vc) 
9(2ei?rc - 1)^ ' 
4e{ercH - H + 2eH^r^) 
3H{2eHrc - I f ' 



(35a) 
(35b) 
(35c) 
(35d) 



From these formulae, it follows that the bulk field fl has 
dimensions of (length) 2. The bulk wave equation (|28p . 
boundary condition (j32[) and (j33[) are the equations we 
must solve. Once we know A and fl the metric pertur- 
bations $ and ^' can be obtained by differentiation. An- 
other quantity of interest is the curvature perturbation 
in uniform density slices, which is given by 



. ^ Ha. 1 ^ 

C = 3' + dq + -A, 

p 3 



(36) 
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coefficient scalar case (matter domination) tensor case (pure tension brane) 



Ai ^epa^^'^k^'^-fi 

A2 

A3 -|eH-'JY'7i - §£^^(4 + 371 -973 +2674) - lePji-^a-S Pa-'^ + ^H{2e - H) 

A4 |e(373 - 272) 

As -ieH-Si 1 

As 1/372 

A7 -2H 

As -iefc*a-^/2^4 

Ag 

Alo 

All = K ^eH-^H'-yi 

A12 = A'lo 



TABLE I: Dimensionless coefficients appearing in tlie canonical wave equations (|25|) for the case of scalar perturbations in the 
late universe and tensor perturbation of a pure tension brane. Even though the An and A12 coefficients do not appear in (I25|) . 
they are crucial for the numeric scheme developed in ijlVI We make use of the 7-factors defined in Eq. (|35|) and the following 
notation: H' = dH /di = —(1 + w)p{l + |/i)~^''^, fi = p + a, and p = p*a~^. 



assuming matter domination. This can be explicitly rep- 
resented in terms of A and Q^,: 

This quantity is interesting because it is expected to be 
conserved on supcrhorizon scales for any metric theory 
of gravity [23j], including the DGP model. Hence, an 
explicit verification that C is constant when k ^ Ha 
provides a useful consistency check of our numerical code. 

We can define a dimensionless canonical bulk field 
by 

1/2 

^ = ^7rT"- (38) 

Substitution of this into (I28|) confirms that -0 satisfies the 
canonical bulk wave equation (I25ap with potential 

y(r,z) = ^^^. (39) 

We can then then replace with -0 in ([5^ and ([55]) and 
make use of (fT8| to obtain the canonical boundary condi- 
tion (j25bp and A equation of motion (|25cp . The explicit 
expressions for the dimensionless A„ coefficients are given 
in Table [J Finally, it is useful to define a dimensionless 
version of as follows: 

Vt = i)^/^ = aZ^r-"^^- (40) 

Generally speaking, it is more useful to work with than 
17 for numeric computations, and later on we will present 
plots of instead of il. 



C. Tensor perturbations 

For tensor perturbations, we restrict ourselves to cases 
where there is no ordinary matter on the brane, 

p = 0, <7 = 3i/re(Hre-e). (41) 

Tensor type perturbations of the bulk geometry are de- 
scribed by perturbations of the form: 

ds^ = -r'^dudv + v'^ {5ij + E^j )dx' dx^ , (42) 

where the 3-tensor Eij is given by 

^-■= E [ 7^^^(-'-y''''4i^)- (43) 

Here, efj is a constant polarization tensor satisfying 

a<,e^^.(k)=5'^e^(k) = fc^(k)=0. (44) 

For brevity, we omit the k subscript and A superscript 
on the mode amplitude E^{u, v) below. Perturbations of 
the bulk Einstein equations yield 

Note that this is equivalent to aiEe'^''') = 0. Per- 
turbation of the junction conditions yields that 

{dyE% = (^E^ + iHE^ + ^^b) . (46) 

To bring this into the canonical form wc introduce the 
new bulk variable 

„3/2 

^ = (47) 
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which satisfies the canonical bulk wave equation with po- 
tential 



(r + z)2 



(48) 



Substitution of ([47]) into the boundary condition ((46l) 
yields the other canonical coefficients A„ in which 
are listed in Table [D Note that since there is no 
brane field A in this case, we set A„ = for n = 
1,2,6,7,8,9,10. 



IV. NUMERIC METHOD 

A. Computational grid 

In this section, we develop a numeric algorithm to solve 
the system of equations (1^51) over a finite region S of the 
(r, z) plane. Illustrations of the shape of E are given in 
Fig. [T] for a few different brane trajectories and choices of 
e. The computational domain has three distinct bound- 
aries: the brane 9Eb, a past null hypersurfaces 
and a future null hypersurface Note how the com- 

putational domain is to the right of 9Sb for the self- 
accelerating case e = -1-1 and to the left for the normal 
case e = — 1. The brane size at the beginning of the sim- 
ulation (where 9Sb and 9E_ intersect, marked i in the 
figure) is a^, while the size at the end (where 9Eb and 
(9I]_i_ intersect) is a/. 

Our numerical algorithm employs an irregular com- 
putational grid as shown in Figs. [T] and [H To define 
this grid, we introduce an arbitrary time parametriza- 
tion along the brane given by the monotonically increas- 
ing function x — x{t). The portion of the brane between 
a — CLi and a = is subdivided into piecewise linear 
segments equally spaced in the new x parameter. That 
is, the change in x over a given segment is 5x = h, where 
h is the overall stepsize parameter of the algorithm. The 
bulk grid is then completed by drawing null lines ema- 
nating from the endpoints of these segments as shown 
in Fig. [TJ As in previous work, the grid involves a num- 
ber of triangular cells adjacent to the brane and diamond 
shaped cells in the bulk. The bulk cells will generally not 
be uniform in size due to our choice of brane partitioning. 

Our actual choice of a; = x{t) is motivated by the desire 
to obtain a quickly-converging algorithm that samples 
the expansion history of the brane sufficiently densely. 
We have experimented with a number of possibilities and 
have found that x{t) — Inrb(i) works best for the self- 
accelerating branch. On the other hand, x{t) = T\^{t) 
seems to give good performance for the normal branch. 
Of course, the choice of x ultimately does not matter 
since all possibilities should give the same results in the 
h ^ limit. 



B. Evolution near the brane 

In order to model the evolution of A and "0 near the 
brane, consider the geometry shown in Fig. [21 Roughly 
speaking, our goal here is to develop an algorithm to 
calculate the values of the fields at the nodes G and H 
given the knowledge of their values at (A,D,E,F). 

If we integrate the bulk wave equation over the trian- 
gular cell EHF, we obtain 



2V'F-'0H-'0E 



A Je J a 



di{XiA + A2S + Ag^Ab + A4< + As^Ab)- (49) 



Here, we have used (j25bp to substitute for the normal 
derivative of ip. We can also integrate (|25cp and (j25dp 
over the brane segment from E to H, we get 



di{XeA + A7S + AsV-b + AgV^I, + AioV'b), 



Ah-Ae= / diE 



(50) 



We now replace the integrals in the exact expressions 
and (|50p with discrete approximations. First, let 
us consider the 2-dimensional integral in (I49p . A simple 
linear approximation to inside EHF yields 

Sp 

d^xV^ = —{VE^pE + VH^H + VF^PF) + 0{Si^). (51) 

A 

Here, 6t is the proper time interval between the nodes E 
and H, which is explicitly 



Si = di= dxT = i/i(T) + 0{h^). 
Je Je 



IE JE 

Here, we have introduced the notation 



{X) = Xh + Xe, 
and defined 

dt dTfj 
dr^ dx 



ixw^x^ 



x^ 



(52) 



(53) 



dx 



7:\Hd 



1 
Hd 



The last equality in ([5^ follows from the trapezoidal ap- 
proximation for one-dimensional integrals. We use this 
same approximation for the other line integrals in ([^5]) 
and ([50)) . after integration by parts to remove the ■0b 
terms and a change of variables from t to x. It is useful 
to simplify the notation by introducing some new coeffi- 



cii = iT. 



cients: 








Cl = 


iTAi, 


C6 = 


^TfAe, 


C2 = 


5TA2, 


C7 = 


iTAy, 


C3 = 


5TA3, 


C8 = 


^TTAg, 


C4 = 


i(A4 - A5), 


cg = 


i(A9 - A'lo 


C5 = 


T^'As, 


ClO = 


T^"^Aio, 



(55) 
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(/5,,a,,o/) = (2.0,0.1,3.0) (p.^,a,,af) = (10.0,0.1,3.0) {a,a,,,af) = (2.0,0.1,5.0) 




bulk spatial coordinate z 
FIG. 1; Typical computational grids used to solve the perturbation equations. 



X = xq + 3h 




FIG. 2: Grid geometry used to derive evolution formulae in 
^TVBl and ^^^C\ The principal brane nodes A, E, H and M 
are separated by a brane time interval 5x = h. We have also 
introduced half-step nodes D, G and L, which are separated 
from the adjacent principal nodes by Sx = h/2. The half- 
step nodes are needed because of the non-local nature of the 
boundary condition in the DGP model. 



that they are known exactly. In terms of these we get 

2V'F - = /i(ciA + csS + csVb + C4Vb> + \\c5ip'\\ + 
^h^{cuf[{Vij) + VfM + 0{h^), (56a) 

||S|| = /l(c6A + CyS -f C8^b + CgV"*) + llcioV'bll 

+0{h^), (56b) 
||A|| = /i(cnS)+0(/i3), (56c) 



where 



d'i/'b dV'b dt 
dx di dx 



(57) 



We can now make use of the following approximations: 



(58a) 



^ _36V.E + 32^D - + 13^H ^ ^ 3)^ (58b) 

to eliminate the ip derivatives in ((56|l to order /i'^. Once 
([58]) is substituted into ([56|) . we have a linear system for 
(tAh, Ah,Sh) in terms of the values of V' at the nodes 
(A,D,E,F), the value of A and 2 at E, and h. Once this 
linear system is solved"^ we know the values of all the 
fields at H (accurate to order h^). Finally, we can use 



V'G = fV'E 



V'D + T^H + T^/'A + 0(/»^), (59) 



to get the value of ■0 at G. Note that the values of A and 
5 at the nodes D and G have not entered the discussion; 
it turns out that it is not necessary to keep track of the 
brane degrees of freedom at non- vertex nodes. 



It is worthwhile noting that all of these Ci coefficients are 
functions of the brane trajectory only, and we assume 



^ Although this is simple to do, the expUcit solution is rather long 
and we omit it from the current discussion. 



8 



C. Evolution in the bulk 

In addition to evolving ?/', A and S near the brane, 
we also need to evolve ip in the bulk. We can take the 
diamond CIJF to be a typical bulk cell. By simply inte- 
grating the bulk wave equation (|25ap over this cell and 
using the divergence theorem, we obtain 



2(V'F+V'i-V'J-V'c)= / d^xVij. (60) 
Jo 

Using a bilinear approximation for the integrand yields 

Vc^c + Vii^i + Vji^F + Vi - i'c)] + 0{h^). (61) 

Here, 5u — 0{h) and Sv ~ 0{h) are the dimensions of 
the cell in null coordinates. Hence, given the knowledge 
of ip on the past nodes F, I and C, we can obtain the 
value at the node J accurate to order h^. 



D. Initial data and computational algorithm 

Having obtained the formulae that tell us how to evolve 
the fields across individual cells, we are now in a position 
to discuss our overall computational strategy. For sim- 
plicity, we will describe how the calculation is carried 
out on the sparse grid shown in Fig. [21 but the method 
is easily generalized to the denser grids shown in Fig. [TJ 

For the numeric solution of a conventional hyperbolic 
problem, it would be sufRcient to specify initial data for 
tp on the nodes (A,B,C,I) in addition to Aa and Ea- 
However, due to the nonlocal boundary condition asso- 
ciated with DGP perturbations, we need to also specify 
"00, ■0E, Ae, and Se initially (see [l^ for a detailed dis- 
cussion of wellposedness and initial conditions for DGP 
perturbations). Once the initial data has been selected, 
the algorithm proceeds as follows: 



1. the diamond evolution formulae (j6ip is then used 
to obtain at the nodes F and J; 



2. the triangle evolution algorithm developed in i jlVBI 
gives 0G, ■'/'H, Ah, and Sh; 

3. Eq. (pT|) is then used to find i/'k; and finally, 

4. the triangle algorithm gives the field values at the 
remaining nodes L and M. 

Obviously, if we have a larger grid than the one shown 
in Fig. [21 steps 2 and 3 need to be iterated a number of 
times. 

For generic grids, the number of diamond cells in the 
grid will scale as while the number of triangle cells 
goes like 1/h. Since the errors involved in the diamond 
and triangle evolution formula are 0{h*) and 0{h^), re- 
spectively, we obtain a final answer that is quadratically 
0{h'^) convergent. 



V. TENSOR PERTURBATIONS ABOUT A DE 
SITTER BRANE 

It is possible to analytically solve the tensor mode per- 
turbation equations given in i jlll CI for the case of pure 
tension branes. The governing bulk wave equation is sep- 
arable in Gaussian- normal coordinates, which allows for 
a Kaluza-Klein mode decomposition of the gravitational 
wave amplitude E. Within the spectrum, one can al- 
ways find a single discrete mode that is normalizablc in 
the bulk. The 4-dimensional mass of this excitation is 
given by [11,111 



'"0 ^ J (Hr^y ' 

lO, 



^ e = -Hi and Hr^ > 2/3, 
e = -1. 



(62) 



In practical terms, the existence of this bound state 
means that for any choice of initial data, we expect the 
late time behaviour of the brane amplitude to be well 
described by a solution of 

+ SHEb + + "^o) -Eb = 0, (63) 

where a = a{t) — e^*. It is easy to verify that this 
implies 

r p m /«"'^'"'% e = +1 and Hr, > 2/3, 
hm Eb{t)(x i (64) 
t-*oo \a", e = -1. 

In Fig. [31 we plot a few typical results for the on brane 
profile of E. All simulations share the property that they 
approach scaling solutions at late time: 



Eb 



(65) 



In all of our simulations for e = — 1, we find that 7 = 0, 
which is consistent with the analytic expectation that 
the normal branch spectrum contains a massless bound 
state. On the other hand, for the self-accelerating brane 
we would expect 
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1 



e = +l. 



(66) 



Also in Fig. [31 we plot this theoretical expectation versus 
our numerical results, and we see excellent agreement. 
Hence, we have confirmed that our code reproduces an- 
alytic results for the gravitational waves in the pure ten- 
sion DGP model. 



VI. SCALAR PERTURBATIONS IN THE 
SELF-ACCELERATING UNIVERSE 

A. Cosmological parameters for DGP late-time 
acceleration 

In this section, we concentrate on the w — a ~ and 
e = -f 1 DGP scenario as a model for the late-time accel- 
erating universe. Examining Tableland Eqns. (|19l23l39p 
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-1 1 2 3 -1 1 2 3 0.5 1 1.5 2 2.5 3 

scale factor log a Hubble constant log Hr^ 



FIG. 3: Numeric solutions for the amplitude of tensor mode perturbations Eh on a pure tension brane (left and center) and 
comparison of the late time power law index derived from analytic and simulation results (right). Note that our late time 
simulation results for e = — 1 are all very similar to the Hr^ = 5 case shown here, and are all consistent with the analytic 7 = 
expectation. 



in detail, we see that the entire evolution of scalar per- 
turbations (modulo initial data) in the canonical formal- 
ism is governed by a single parameter /5*. Recall that 
we defined p* to be proportional to the matter density 
when the mode being modeled crosses the Hubble hori- 
zon. However, we note that in any model which asymp- 
totes to de Sitter space in the future, any mode that 
enters the horizon must exit it as well. Indeed, by solv- 
ing the equation k/ Ha = 1, we find that there are always 
two horizon crossing epochs: 



1 and 



3 + V9 + 12/5, 



(67) 



If /5, > 6, the second solution is always greater than 
the first; i.e., the perturbation enters the horizon at a = 
1 and leaves later. However, if = 6 both solutions 
coincide and horizon entry and exit occur at the same 
time. It is easy to confirm that the physical wavelength of 
the /5, = 6 mode coincides with the Hubble length at the 
moment when the brane switches from the decelerating 
to the accelerating phase; i.e., when a = 0. 

We can go further by incorporating data from cosmo- 
logical observations. By examining probes of the expan- 
sion history [l^, it has been found that 



1 



0.15 ±0.02, 



(68) 



at 95% confidence.'* Once the value of Vlr^ is fixed, we 
can use it to find the value of a, for any given mode in 
terms of : 



1 3(1 -2f7 



1/2. 



(69) 



* Actually, in [Tsl l the constraint Q,n 
which I I68I I can be easily derived. 



0.23 it 0.04 was given, from 



This in turn allows us to determine the comoving 
wavenumber from 



= 2nl{^ka^HQ, iJo^i = 2998/1-1 Mpc. (70) 



Hence, we can find k in terms of p*. In order to facilitate 
the comparison of our results with observations and the 
literature, it is more convenient to invert the procedure 
to obtain p.^, in terms of k. However, this only works if k 
is above a critical value 



fee = 2nll^{l - 2nll^f'^Ho = 0.0003 /iMpc 



(71) 



The last equality follows from taking the best fit value 
for r^r^, which is what we do from now on. Any modes 
with k < kc have physical wavelengths larger than the 
Hubble length throughout the cold dark matter and late 
time acceleration epochs, and therefore are not modeled 
by our code. Finally, it is sometimes useful to have an 
explicit expression for the gravitational master variable 
rib in terms of the canonical field iph- This is: 



-3/2 



(72) 



B. Typical waveforms 

In Fig. m we plot the results of our simulations 
for several values of k greater than the critical value 
0.0003 /iMpc^i. For all plots of scalar perturbations in 
this paper, we select the bulk field to be zero and the 
brane field non-zero initially. We have also simulated 
several different choices of initial data, such as the bulk 
field being constant along the initial null hypersurface, 
and have found that the simulation results remain the 
same as long as the initial time is early enough. This is 
analogous to what happens in the RS case 11 1 . 
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From Fig. [4] it can be seen that we recover ordinary 
4-dimensional GR at very early times. In particular, for 
a ^ a* and all values of fc > /cc we see that: 

• the metric perturbations are conserved and have 
the opposite sign, $ w — 

• the density perturbation is proportional to the scale 
factor, A oc a; and, 

• the bulk master variable scales as fib oc a*^ on the 
brane. 

In addition, we have checked that the 4-dimcnsional Pois- 
son equation is satisfied before horizon entry: 



In addition, the following relations are predicted to hold: 



1 



(73) 



In other words, we have explicitly confirmed that DGP 
perturbations behave as in GR on superhorizon scales 
before horizon crossing. 

Finally, in Fig. [5] we plot the simulation results for the 
behaviour of the curvature perturbation C as given by 
([37|) for a few different large scales. As can clearly be 
seen, ( is conserved for both early and late times when 
the physical wavelengths of the modes are much larger 
than the horizon size. This is to be expected for any 
metric theory of gravity [23| , and hence provides a good 
consistency check of our code. 



$ = 



? 9 



1 + — A. 



(78a) 
(78b) 



In Fig. [6l we compare simulation results versus the 
QS approximation for the linear growth factor g{a) = 
A{a)/a, and the alternate gravitational potentials <E>± — 
i($ ± ^P). We see that the simulation results are consis- 
tent with the QS approximation for k > 10^^ hMpc'~^ . 

We further quantify the performance of the QS approx- 
imation in Fig. [71 There, we show the relative error in 
the QS approximation as a function of the scale. This is 
defined by 

QS prediction — simulation result 



simulation result 



rel. error = — ^ 100%. 

(79) 

We see that the relative error in the QS prediction for 
A is fairly low (< 4%) on all scales. Conversely, the QS 
values of $± become reliable only for k > 0.01 /iMpc~^, 
with errors of less than ~ 5%. 



D. Superhorizon behaviour in the asymptotic 
future 



The quasi-static approximation and subhorizon 
behaviour 



In [81], a 'quasi-static' (QS) approximation was devel- 
oped to describe the behaviour of DGP perturbations 
whilst well inside the cosmological horizon and with phys- 
ical wavelengths much less than the crossover scale: 



k ^ Ha, a <C krc 



(74) 



These conditions wih hold for modes with k ^ kc and 

1 /2 

k > 2nri Ho (up to some redshift); or, equivalently, if 



k > 10"'*/iMpc" 



(75) 



In this section, we compare the QS approximation to our 
simulations to determine just how large k must be for it 
to be valid. 

In the QS approximation, one neglects the time deriva- 
tives of 17 compared to the spatial gradients. This allows 
one to solve the bulk wave equation (j^B]) , and hence close 
the system ([32]) and ([33l) on the brane. This leads to the 
following ordinary differential equation for A: 



A + 2HA 



1 



1 



where 



P=l- 2eHrc 1 



1 

3(3 

H ' 

3IP 



(76) 



(77) 



In this subsection, we attempt to explain/predict the 
very late time behaviour of our simulations by demon- 
strating that there exists a bound state of the bulk field 
in the asymptotic future of the evolution. As we have al- 
ready seen for the case of tensor perturbations in ifVl such 
bound states tend to dominate the late time behaviour 
of the model, irrespective of initial data. 

In the asymptotic future, the brane geometry ap- 
proaches that of de Sitter space with H — I/tc- In the 
de Sitter regime, the bulk wave equation ([28)1 becomes 



3 dVl 
Tc dt 



2 



y\dn 
TcJ dy 



ril, (80) 



with a — e*^^". This equation is solvable via the separa- 
tion of variables Q{t,y) — T\(t)ujx{y), where 



dy 



3 dTx ^ 

Tc dt 
1 dhJ\ 



a2 r? ' ^' 



{y + r^Y dy 



A 



(81a) 
(81b) 



Here, A is a dimensionless separation constant. The so- 
lution for L0\ is 

uJx = a+(l + ^\ +a-(l + ^\ , (82) 
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FIG. 4: The results of our simulations on the brane for several choices of k. We have normalized the value of $ to be unity 
at early times. Also note that the lower left panel shows the dimensionless bulk master variable f2b, as defined in Eq. (140 |l . 
divided by a^. All simulations are performed with p* > 6, which means that all modes enter the horizon at a = 1, or when 
a = a,. 
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FIG. 5: Behaviour of the ( curvature perturbation on large 
scales for the self-accelerating branch (we have normalized 
= 1 at early times) . Note that the curvature perturbation is 
conserved when the modes are superhorizon; i.e., at both early 
and late times. This is to be expected for any conservative 
theory of gravity, such as the DGP model. 



where iy± = | (^1 ± y/l + |Aj . Assuming that A is real, 
we need to set a-)_ = to ensure that lox is normalisable 



under the Sturm-Louisville inner product; i.e., that 



dy 



(y + r-c)4' 



(83) 



is finite, which means we have a true bound state. 

Now, if we put H = l/r^ in the boundary condition 
([5^ we obtain 



r.1 



1 



A 



I'q \ U" / p / K" 

(84^ 

Let us now assume that for very late times (fcrc ^ a) the 
following conditions hold: 



fc2 



A 



< min ( |f^b|, r^. 



\n 



(85) 



i.e., the A term is negligible on the righthand side of ((84| . 
Under these assumptions, which we still need to justify, 
([80| and ((84|) form a closed system for fi. Putting our 
mode solution fl{t,y) — Tx{t)uj\{y) into the boundary 
condition with = and neglecting A, we obtain 



A 



(86) 



Putting A = — 2 into the temporal equation (|81ap and 
solving for T\, we find that 



Tx{t) ^bia^{t) + b2a{t), 



(87) 
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FIG. 6: Linear growth factor and alternate gravitational po- 
tentials "1>± from simulations and the QS approximation in 
the self-accelerating branch. In the top panel, we normalize 
g{a) to unity at early times, in the lower two panels we nor- 
malize "I>_ to unity as a — > 0. For comparison, we also show 
the relevant results for the concordance ACDM model with 
Vim = 0.26 and = 0.74. 



where hi and 62 are constants. Of course, the 61 solution 
will eventually dominate, which leads to the following 
asymptotic bound state: 



nn 1 



where fto is a constant. However, before we assume that 
the late time behaviour of the system is indeed described 
by this bound state, we must verify that the assumptions 
([55)1 under which it was derived are valid. To do so, we 
note that when H = 1/rc, the A equation of motion 
reduces to 



A-^ — A- 



kIpA 



3a5 



r!b. 



(89) 




wavcnumbcr log k [h Mpc 



FIG. 7: The relative error in the QS approximation, as defined 
in Eq. (|79|l . for various quantities evaluated in the present 
epoch and assuming the self-accelerating branch. For k > 
0.01 hMpc~^ the errors are less than 5%. 



relevant, so we just quote the late time behaviour 



Ait) ^ Ao, 



(90) 



where Aq is a constant. With the solutions ([55]) and (|90p . 
we see that the assumptions ([55)1 are indeed satisfied at 
sufficiently late time. Hence we have succeeded in finding 
an asymptotic bound state that is expected to dominate 
the system's behaviour at late time. Finally, note that 
we can use these asymptotic solutions for 57 and A with 
(l34t to obtain 



f2oa(i) 



(91) 



Making use of ([58)) . this equation can be solved exactly. 
However, the full solution is complicated and not really 



i.e., $ « 5" at late time. We have verified that the asymp- 
totic solutions (|88l90l9ip are realized in our simulations 
at late times. 

Before moving on, we would like to remark on the ap- 
parent instability of the self-accelerating DGP model as 
indicated by the divergence of 17, $ and 5* in the asymp- 
totic future. This unstable mode corresponds to the ra- 
dion, which is a physical degree of freedom in the self- 
accelerating branch despite the fact that we have only 
one brane [12j |. It is well known that in this case the ra- 
dion has a negative mass squared — —47?^ and thus it 
is unstable [2y]. However, as was shown in Ref. [13], this 
is not a true gravitational instability on the brane as it is 
possible to find a gauge in which all metric perturbations 
remain finite. 
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VII. SCALAR PERTURBATIONS IN THE DGP 
NORMAL BRANCH 

A. Cosmological parameters for ADGP 

We now turn our attention to the behaviour of density 
perturbations in the normal branch of the DGP model. 
Unlike the e = +1 case, this branch does not naturally 
have a late time accelerating phase. So in order to be 
made consistent with observations, we must allow for 
the brane to have a nonzero tension that acts as an ef- 
fective 4-dimensional cosmological constant (we call this 
the ADGP model). Assuming that the matter sector is 
CDM-dominated, the Friedmann equation for this sce- 
nario follows from the general form with e = — 1 and 
w — 0. The background d yna mics has been compared 
to observations of H{z) in [l^|, who finds the following 
parameter values: 



0.23 ±0.04, nr 



< 0.05, (92) 



at 95% confidence. Here, po is the present day CDM 
density. Note that the observationally preferred value of 
^Ir^ is zero. Since the DGP model goes over to GR in 
this limit, this implies that ACDM gives a better fit to 
the data than ADGP. In what follows, we will always 
assume the best fit value of 0.23 for ilm and treat ^Ir^ as 
an adjustable parameter that must be smaller than 0.05 
to yield a reahstic model. 

Once rim and have been selected, it is straight- 
forward to obtain the value of the dimensionless brane 
tension: 

3(1 - + 2nl.i^) 



2 2 



(93) 



which can be re-written in terms of a new density param- 
eter rirr- 



p 
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FIG. 8: Linear growth factor and alternate gravitational po- 
tentials "1>± from simulations and the QS approximation in the 
normal branch with 57,.^ = 0.05 and 0,-^ = 0.23. As in Fig. |6] 
we normalize 17(a) and <&_ to unity at early times. For com- 
parison, we also show the relevant results for the concordance 
ACDM model with 0.,-^ = 0.26 and = 0.74. 



0„ 



3i/2 



1-f^n 



(94) 



One can also find a* in terms of the observational param 
eters and 

Q 1 sr^ni 



(95) 



As before, these two formulae can be used in (|70p to de- 
termine k as an explicit function of , or as an implicit 
function of k. Following the latter approach means that 
when we select k along with fim and fi^^ , the evolution 
of perturbations is completely specified up to the choice 
of initial data. 



B. Simulation results and comparison to the QS 
approximation 

In Fig. [51 we compare the results of our simulations to 
the QS approximation and AGDM in the case = 0.05. 



As in WIBl wc find that the simulation results are fairly 
insensitive to initial conditions provided that the initial 
data surface is set far enough into the past; here, all 
plots have been generated assuming f2 = initially. In 
contrast to the self-accelerating case, we find that the 
linear growth factor and $_ potential are generally larger 
than in the AGDM case. The general trend is for <!>_ to 
become larger on small scales. We also notice that the 
QS approximation seems to provide a very good match 
to the simulation results for A on all scales. 

In Fig. [51 we show the effect of changing the fi^c param- 
eter on the simulation results for <&_ . For any given scale, 
we see that the £7^^ limit approaches the ACDM 
prediction. Also, we note that the simulation results are 
closer to the ACDM case for smaller values of k] i.e., the 
most pronounced deviations from GR are observed on 
the smallest scales simulated. 

Finally, in Fig. [TU] we quantify the error in the QS 
prediction for the value of various quantities at z = 
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FIG. 9; The TSW potential' $_ as a function of the scale factor for various scales and choices of Qr^ for the normal branch. 
In all cases we have taken = 0.23. The ACDM curves are included for purposes of comparison. 



as a function of the scale. As in the self-accelerating 
case, we see that the QS approximation provides rea- 
sonably accurate results (with errors < 5%) on scales 
k > 0.01 /iMpc"^ 



C. Superhorizon behaviour in the asymptotic 
future 

In the asymptotic future, the brane geometry ap- 
proaches that of de Sitter spacetime with H determined 
by cr 7^ 0. Unlike in the self- accelerating branch, there 
appears a horizon at y = 1/ H . An analysis similar to the 
one presented in WI Dl for the self-accelera ting branch 
shows that there is no solution with a real A [23] ■ There- 
fore, there is no bound state solution in the asymptotic 
de Sitter spacetime and f2 is a superposition of massive 
Kaluza-Klein modes that oscillate in time. It is worth 
noting that the dynamical scaling ansatz implicitly as- 
sumes the existence of a O bound state, and will hence 
fail in the asymptotic de Sitter future of the DGP normal 
branch. 



VIII. CONCLUSIONS 

In this paper, we have presented numeric solutions 
for cosmological perturbations in the DGP braneworld 
model both in the self-accelerating and the normal 
branches. We extended the algorithm developed for the 
Randall-Sundrum (RS) model to handle the nonlocal 
boundary conditions characteristic of the DGP model. 
The numerical code was tested for tensor perturbations 
and the agreement with the analytic solutions was found 
to be excellent. 

We confirmed that on small scales k > 0.01ft. Mpc~^, 
the quasi-static (QS) approximation reliably predicts the 
evolution of perturbations with relative errors less than 



around 5% at z = 0. Our results are quite insensitive 
to the initial conditions as long as we start our simula- 
tions early enough. On larger scales, the potential 
which determines the integrated Sach- Wolfe (ISW) effect, 
shows more suppression than the QS prediction. We find 
that our numerical solutions agree well with the dynam- 
ical scaling (DS) solution both in the self-accelerating 
and normal branches, except in the asymptotic de Sitter 
phase of the normal branch where the dynamical scaling 
solution fails to exist [2^ . 

Our numeric solutions provide the basis for studying 
observational signatures of the model, especially in the 
normal branch where the influence of the extra dimension 
on the evolution of large scale structure has not yet been 
explored. 
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APPENDIX A: BOUNDARY CONDITION AND A 
EQUATION OF MOTION FOR SCALAR 
PERTURBATIONS 

To derive the boundary condition satisfied by the bulk 
master variable f2 and the second order equation of mo- 
tion for the density perturbation A for scalar perturba- 
tions, we begin with the linearized version of the effective 
Einstein equations d?]): 

SG^"^} = (2Kln)Hn^, - S£^,. (Al) 
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FIG. 10: The relative error in the QS approximation, as 
defined in Eq. (|79p . for various quantities evaluated in the 
present epoch and assuming the normal branch. For k > 
0.01 /i Mpc~^ the errors are less than 5%. Also notice how 
the QS approximation is generally more accurate for smaller 



In this expression, SG\^} and 511^,^ can be obtained by us- 
ing the perturbed brane metric (|29p and perturbed stress 
energy tensor l|30|) . We can parameterize the perturba- 
tions of the bulk Weyl "fluid" Ef^^, as follows: 



aSq, 



'£,i 



\5pe 5\ + 5t:, 



(A2) 



Here, Suf^ — Stt^^j — ^Stt^j^'^ 5ij , a comma denotes partial 
differentiation and indices are raised and lowered with the 
flat 3-metric. Using this in the (0, i) component of the 
perturbed effective Einstein equations (|Aip . we obtain 
the following equation: 



iJ* - $ = ^ 



2H7\ 



2Hi\ 



pV 
k 



2Hrce 



(A3) 



where we have made use of the fact that the background 
brane matter distribution is CDM plus a possible effec- 
tive cosmological constant induced by the brane tension. 
Combining this with the (0, 0) component of (|A1|) yields 
the Poisson equation: 



2eHrc - 1 



pA 



5p£ - 3HSq£ \ 



2eHrc 



(A4) 



In these formulae, the gauge invariant density pertur- 
bation A is defined in (|3ip . while the invariant velocity 
perturbation is given by V — —kSq/pin the longitudinal 
gauge. 

Additional relations can be obtained by noting that 
5(V"r„0) = 0. (A5) 
The spatial and temporal components of this yield that 



V + HV = (A6a) 
a 

A = -^(^1-^H^V-3{^^H^), (A6b) 

respectively. Combining (|A3p with (jA6p yields a second 
order differential equation for A: 



P 3 
A + 2HA^ — T* + 7:^ 
2 



3HF, 



with 



F = 



2Hrce - 1 ' 



(A7) 



(A8) 



Now, one important feature that distinguishes 
braneworld cosmological fluctuations from the GR case 
is that in addition to perturbations of the geometry and 
the matter, we must also consider perturbations of the 
brane's position. That is, in the Gaussian normal co- 
ordinates system the brane is located a.t y — before 
perturbation and at y = ^ after perturbation, where ^ 
is the scalar brane bending degree of freedom. It is use- 
ful to parameterize the perturbed geometry of the y = 
hypersurface (i.e., the brane's unperturbed position) by 

dsl^o = -(1 + '2A)dt^ + a^il + 2n)dijdx'dx^ . (A9) 

The metric potentials at the unperturbed brane position 
{A, TV) are then related to the metric potentials at the 
perturbed position (^f, <I>) by 




U, $ = 7^ - eH£,. (AlO) 



Deffayet [2^ has shown that the brane bending scalar is 
simply given by 



-re($ + *). 



(All) 



In addition, he demonstrated that it is possible to express 
A and TZ in terms of the bulk master variable f2: 



6a 



ba 



H 



2fc2 



3e — - i? {dyn)b + -^rib - 317b + eij^t 



H 



3ei/((9yrJ)b + ^r^b - SHili, 



(A12) 
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It can also be shown [221 that the Weyl fluid perturba- 
tions are also directly given by O: 

= '^l^is^-^iHn-n). (A13) 

Now, these formulae in conjunction with the wave equa- 
tion (I^S)) can be used to re-write the Poisson equation 
(PI)) as: 

ju2 1.2 

2r,KlpA = 4rc^<i>+— e -[idyQ)b~eHn^]. (A14) 

a"^ a-^ 



Then, one can use (jA10p - (jA12p in this equation to obtain 
the boundary condition (|32p . One can then use ([221) with 
(|A10p - (|Al2|l to get $ and ^> in terms of fib and A; i.e., 
Eqs. ([Ml). Finally, substituting ^ and (|XT3)l into ([Tzl) 
gives the final form of the A equation of motion ([33l) . 
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